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We use the recently proposed directed-path algorithm to study the chiral limit of strongly coupled 
lattice QCD with staggered quarks at finite temperatures. The new algorithm allows us to com- 
pute the chiral susceptibility and the pion decay constant accurately on large lattices for massless 
quarks. In the low temperature phase we find clear evidence for the singularities predicted by chiral 
perturbation theory. We also show convincingly that the chiral phase transition is of second order 
and belongs to the 0(2) universality class. 
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INTRODUCTION 

One of the outstanding problems in lattice QCD is 
to compute physical quantities reliably when the quarks 
have a small mass. All conventional algorithms suf- 
fer from critical slowing down as the quark masses are 
lowered. Today, most calculations use quarks that are 
heavy and the results are then extrapolated to the chiral 
limit using chiral perturbation theory (ChPT). Whether 
ChPT is applicable to the data in the range of quark 
masses that are accessible today, is highly debatable p|. 
The evidence for chiral singularities predicted by ChPT 
is weak. The chiral singularities are logarithmic in four 
dimensions and power-like in three dimensions. Even 
the power-like singularities, which are well known in spin 
models 2], have yet to be seen in lattice QCD calcula- 
tions. 

It has been predicted that QCD with two massless fla- 
vors of quarks will undergo a finite temperature second 
order chiral phase transition in the 0(4) universality class 
1^ ,41 • Although there is clear evidence for such a phase 
transition from lattice calculations, the above difficulties 
with chiral extrapolations also affect our ability to estab- 
lish the universal properties of the phase transition. In 
particular no precision calculations for the critical expo- 
nents exist that match with expectations and rule out 
other universality classes with similar exponents. For ex- 
ample, lattice QCD with staggered quarks contains an ex- 
act 0{2) chiral symmetry which is a subgroup of the full 
chiral symmetry. In a two flavor theory, this symmetry is 
expected to be dynamically enhanced in the continuum 
limit to the 0(4) symmetry. Ideally it should be possible 
to show that on coarse lattices the chiral phase transition 
belongs to the 0(2) universality class and as the contin- 
uum limit is reached the universality class must change to 
0(4). At present results from lattice simulations do not 
appear to match with either 0(2) or 0(4) universality 
class 0,01. Although, results with two flavors of Wilson 
quarks appear to show consistency with 0(4) critical be- 
havior I3, these lattice fermions do not even possess the 
relevant chiral symmetry. The presence of a parity- flavor 
broken phase nearby complicates matters further. More 
work is necessary before one can understand the results 
of 0. 



Given these difficulties it is useful to find at least some 
point in the phase diagram of lattice QCD where preci- 
sion calculations with massless quarks are possible. In 
this article we consider the strong coupling limit of lat- 
tice QCD with staggered quarks. Although, this limit 
has the worst lattice artifacts, it is a good toy model that 
shares some qualitative features of QCD. The quarks are 
confined and the system is known to break the remnant 
0(2) chiral symmetry. Computationally the theory sim- 
plifies enormously since all the gauge field integrals can 
be performed exactly. However, the remaining dynamics 
of quarks is still non-trivial and leads to an interesting 
theory. The phase diagram in the temperature versus 
baryon density plane is also expected to be interesting 
8] . The strong coupling limit was first studied with mean 
field methods 0, 0, Later numerical simulations 

were proposed to understand the theory from first prin- 
ciples Unfortunately, all Monte Carlo algorithms 
developed so far have suffered from critical slowing down 
near the chiral limit. For this reason calculations were 
performed away from the chiral limit which limited their 
precision in determining chiral quantities. The dream 
to accurately solve lattice QCD, even in this simplified 
limit, remains unfulfilled. 

Over the last decade a revolution has occurred in the 
field of Monte Carlo algorithms. A variety of classical 
and quantum lattice models can now be solved accurately 
using cluster algorithms that can beat critical slowing 
down very efficiently. A recent review of the progress 
can be found in 14). Recently, an extension of these ideas 
has led to the discovery of the directed-path algorithm for 
studying the chiral limit of strongly coupled lattice gauge 
theories • For the first time this algorithm allows us to 
precisely compute quantities in the chiral limit without 
further approximations. 

In this article we apply the new algorithm to study 
the finite temperature chiral phase transition in strongly 
coupled lattice QCD with massless staggered quarks. We 
focus on the physics of pions and the universal properties 
near the phase transition at zero baryon density. We use 
U{3) gauge fields instead of SU (3) in order to avoid ineffi- 
ciencies in the algorithm due to the existence of baryonic 
loops in SU (3). The distinction between U (3) and SU (3) 
should not be important for our study since the baryons 
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are expected to have a mass close to the cutoff. We show 
that some of the predictions from three dimensional 0(2) 
ChPT are borne out in the low temperature phase of the 
model and establish with high precision that the chiral 
phase transition belongs to the 0{2) universality class. 



THE MODEL AND OBSERVABLES 

The partition function of the model we study in this 
article is given by 

Z(r, m) = J [dU] [d^dip] exp {-S[U, ip, {p]) , (1) 

where [dU] is the Haar measure over U{3) matrices and 
[dipdip] specify Grassmann integration. At strong cou- 
plings, the Euclidean space action S'[?7, t/), ■0] is given by 
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(2) 

where x refers to the lattice site on a periodic four dimen- 
sional hyper-cubic lattice of size L along the three spatial 
directions and size Lt along the euclidean time direction, 
11= 1,2,3,4 refers to the four directions, Ux^/j, is a 3 x 3 
unitary matrix associated with the bond connecting the 
site X with the neighboring site x + p, and represents the 
gauge field, ipx is a 3-component column vector and ip^ 
is an 3 component row vector made up of Grassmann 
variables and represents the staggered quark field at the 
site X. We will assume that the gauge links satisfy pe- 
riodic boundary conditions while the quark fields satisfy 
either periodic or anti-periodic boundary conditions. The 
factors r]x^fj, are the well known staggered fcrmion phase 
factors. We will choose them to have the property that 



X,fJ, 



1,^ = 1,2,3 (spatial directions) and rj^ 4 = T 
(temporal direction), where the real parameter T acts 
like a temperature. By working on asymmetric lattices 
with Lt « L and allowing T to vary continuously, one 
can study finite temperature phase transitions in strong 
couphng QCD 0. 

The partition function given in eq.(^ can be rewrit- 
ten as a partition function for a monomer-dimer system, 
which is given by 

(3-6.,,)! J. 



Z(T,m) 



[n,b] x.fj. 



(3) 

and is discussed in detail in |l2lll5l |. Here Ux = 0, 1, 2, 3 
refers to the number of monomers on the site x, bx,^ = 
0,1,2,3 represents the number of dimers on the bond 
connecting x and x -\- fi, m is the monomer weight, 
Zx,iJ. = Vx are the dimer weights. Note that while 
spatial dimers carry a weight 1/4, temporal dimers carry 
a weight T/4. The sum is over all monomer-dimer con- 
figurations [n, b] which are constrained such that at each 
site, Ux + J2n[bx,^L + bx-fL.,fj] = 3. 



The model described by Z{T,m) is known to have an 
exact 0(2) chiral symmetry when m = 0. This symmetry 
is broken at low temperatures but gets restored at high 
temperatures due to a finite temperature chiral phase 
transition. In order to study the chiral physics near this 
transition we focus on three observables as defined below: 



(i) The chiral condensate 



L-^ Z am 



(ii) The chiral susceptibility 
1 1 92 



X 



(4) 



(5) 



and 



(ii) The hclicity modulus 



where J.,, ^ 
sites and 



(Txibx.fj^ - N/8), with ax 
= — 1 on odd sites. 



1 on even 



When m = the current Jx,^ is the conserved current 
associated with the 0(2) chiral symmetry. Further, as 
discussed in 0|, it can be shown that = \imL-,aD Y, 
where F is the pion decay constant. 



UNIVERSAL PREDICTIONS 

Let us now briefly review the universal predictions rel- 
evant to our study. The predictions from chiral pertur- 
bation theory for 0{N) models have been discussed in 
|l6l | . In particular the finite size scaling formula for x at 
TO = is given by is given by 



X 



— V2r3 



1 



l + f3,{N-l)— + ^ + 



(7) 



where iV = 2 in our case, ^ Y, f3i = 0.226..., E = 
lim,„_to liiiiL-*oo (</>)) and a is a constant dependent on 
other low energy constants. 

If the chiral phase transition is second order then the 
value of S, computed in the low temperature phase using 
eq. |7J), is a function of the temperature and satisfies the 
relation 



S(T) = A(T, -T)'3,T<r„ 



(8) 



close to the critical temperature Tc- In the critical region 
we also expect x to satisfy the scaling relation 



X 



(9) 



where t = (T/Tc — l) and g{x) is an analytic function of x 
with properties g{x) = go + gix + ... for small x, g{x) — > 
|a;p^ for x —00 and g{x) a;-''(2-')) for x 00. 
The 0(2) universality predicts predicts f3 = 0.3485(2), 
v = 0.6715(3) and r] = 0.0380(4) [l^. 
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RESULTS 

We have done extensive simulations on lattices ranging 
from L = 8 to L = 192 with fixed if = 4 at various values 
of T. We will first verify that our results satisfy eq. (0 
in the chirally broken phase and then show that our data 
in the entire critical region is consistently described by 
eqs.(|Sl) and © with 0(2) critical exponents. 
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FIG. 1: Plots of X vs. L (top left), vs. L (bottom 

left),x - E^LVSfl + I3i/{F'^L)] vs. L (bottom right), and Y 
vs. L (top right) at T = 7.42. The solid lines are fit to the 
data as discussed in the text. 

In order to verify eq. lO, let us focus on T = 7.42, a 
temperature which is below but close to T^. We believe 
that the closeness to the second order transition point will 
reduce lattice artifacts and help make connections with 
ChPT easier. In the top-left graph of figure ^ plot 
our results for % as a function of L. The availability of 
results at large values of L with small error bars allows us 
to fit our results to the first three terms in eq. {T)) reliably. 
Using the values of x for L > 24 we find S = 1.079(2), 
F = 0.181(4) and a = 114(4) with a xVd.o.f of 0.73. In 
order to convince the readers that the error bars in our 
data are small enough to be sensitive to the higher order 
finite size effects, we also plot % — S-^L'^/2 versus L in the 
bottom-left graph and x - 2^^3/2(1 + (3i/{F'^L)) versus 
L in the bottom-right graph of figure ^ As the graphs 
indicate, the evidence for the and L dependence of x is 
indeed quite strong. Since = MmL^oa Y : we can also 
compute F independently using the helicity modulus. In 
the top right plot of figuren]we plot F as a function of L. 
A fit indicates that the data for L > 64 is consistent with 
F'^ = 0.0327(2) (sohd line in the figure). This gives F = 
0.181(1) in excellent agreement with the value obtained 
above. Thus, we find that the L dependence of x at m = 
and T = 7.42 is indeed consistent with expectations 
from ChPT. 

Let us now turn our attention to the entire region near 
Tc. Clearly, the fits based on ChPT (eq.Q) wiU be re- 



liable only for L >> ^, where ^ represents the diverging 
correlation length at Tc- In three dimensions ^ ~ l/F"^. 
Since most of our computations involve lattice sizes with 
L < 128 we assume l/F"^ must be less than 64 for our fits 
to be reliable. With this criterion we find that the above 
fitting procedure is reliable only when T < 7.45. The 
fits of our data in this range of temperatures are given 
in tabled The values at T = 7.42 have been discussed 






E 


F 


a 


xVd.o.f 


7.30 


1.599(5) 


0.25(2) 


6(12) 


0.95 


7.33 


1.502(6) 


0.26(2) 


34(8) 


1.5 


7.40 


1.198(11) 


0.18(2) 


62(11) 


0.1 


7.42 


1.079(2) 


0.181(4) 


114(4) 


0.73 


7.43 


1.014(3) 


0.174(6) 


142(6) 


1.7 


7.44 


0.931(3) 


0.149(7) 


168(18) 


0.15 


7.45 


0.827(9) 


0.127(10) 


242(25) 


0.75 



TABLE L Fits of x vs. L to eq. Q. 



above. We can now verify if the values of S given in ta- 
ble ^satisfy eq. lO. However, we find that using Tc as a 
free parameter in the fits is not ideal since the fits do not 
converge. Thus, we first determine Tc using a different 
method. 

Although, ChPT is unreliable on small lattices close 
to the phase transition, it can still help in roughly es- 
timating Tc. We find that 7.476 < Tc < 7.478. In 
this region we have data at four different temperatures: 
7.4765, 7.477, 7.4775, 7.478. Since for these values of the 
temperature we expect {T/Tc — 1)L^/'^ to be small, the 
0{2) scaling suggests that in this temperature range x 
must satisfy 

X^goL''''+giiT/Tc - 1)L3-45i, (10) 

obtained from eq. © after substituting the 0{2) critical 
exponents. Indeed, a joint fit to all the available data in- 
volving 22 data points yields go — 14.69(2), gi — —9.7(3) 
and Tc = 7.47739(3) with a xVd-O-f of 0.74. This joint 
fit is shown in figure [21 Interestingly, if we use the 3d 
Ising critical exponents in this fit instead of the 0(2) ex- 
ponents, we obtain Tc = 7.47734(3) with a x^/d.o.f of 
about 1. Thus, clearly we cannot distinguish between 
Ising and 0(2) exponents with this approach. However, 
we believe we can determine Tc quite accurately. 

Having obtained Tc we can now verify if the values 
of S given in Table D for various values of T, satisfy 
eq. (0. In figure 13 we plot our results for S as a func- 
tion of T. We find that our data fits well to eq. © 
with Tc = 7.47739 fixed. We obtain A = 2.92(2) and 
f3 = 0.348(2) with a xVd-O.f of 0.53. Changing Tc to 
7.47734 has negligible effect on this result. The value of 
(3 is in excellent agreement with 0.3485(2), the exponent 
obtained in the three dimensional XY spin model 0| 
and differs significantly with the three dimensional Ising 



4 



2x10 



_ 3.0x10 



2.9x10 



1x10 




^0 



40 



60 



80 
L 



o 


T = 


7.4765 


o 


T = 


7.4770 


o 


T = 


7.4775 


o 


T = 


7.4780 



100 



120 



200 



150 



100 



50 



-20 



-10 



-250 -200 



-150 

(T-TJ L' 



100 -50 

(1/v) 



FIG. 2: Plot of X vs. L for four different values of T near Tc. 
The solid lines represent the function given in eq. HIO^ with 
go = 14.69 and gi — —9.7. The inset shows how the fitting 
function fits the data at L = 48. 




FIG. 3: Plot of S vs. T at m = and (0) vs. m at T = Tc 
(Inset). The solid lines represent fits discussed in the text. 

exponent 0.32648(18) [3. In order to test the scaline; 
relation eq.® we plot x/L^"'' against (T - Tc)L^^'' in 
figure^ In order to be fair we do not include the data 
at T = 7.4765, 7.477, 7.4775 and 7.478, since these were 
already used to fit to this relation while determining Tc- 
The remaining data points (a total of 77) fall consistently 
on a single function as shown in figure 01 

So far we had focused on the physics at m = 0. 
The algorithm discussed in can easily be extended 
to include a quark mass. At T = Tc one expects 
limL^oo(0) = Bm^/\ where 6 = 4.780(2) in the case 
of 0(2) universality [l2|. Since this is another indepen- 
dent test of the critical behavior and our prediction of Tc, 
in figure 13 (inset) we also show limi^oo(0) as a function 
of TO at r = Tc = 7.47739. The data fits well to the 



FIG. 4: Plot of xZ-t'^"" vs. (T - Tc)L^^'' . The plot confirms 
the scaling prediction of eq.Q. 



expected form with B = 5.9(1) and S = 4.97(10) with a 
X^/d.o.f of 0.34. On the other hand fixing S = 4.78 in the 
fit gives B — 6.18(2) and increases the x^/d.o.f to 0.88. 



CONCLUSIONS 

The above results show convincingly that the finite 
temperature chiral phase transition, in strongly coupled 
lattice QCD with staggered quarks, is governed by 0{2) 
universality class. Further, close to the transition in the 
low temperature phase the singularities of 0(2) ChPT 
arising due to Goldstone pions is observable. 

There are many directions to extend the current work. 
In the low temperature phase one can verify the striking 
predictions of ChPT for the quark mass dependence of 
the chiral condensate, the pion mass and the pion de- 
cay constant. Using the pion decay constant as a phys- 
ical scale it would also be interesting to understand the 
width of the critical region in physical units. This may 
shed some light on how difficult it would be to study 
the universal properties of the chiral transition at weaker 
couplings. The physics of the chiral limit of SU {N) gauge 
theories at strong couplings in the presence of a baryon 
chemical potential is quite rich and has not yet been re- 
liably explored from first principles. In particular when 
N is even it is well known that the inclusion of a baryon 
chemical potential does not lead to sign problems. Al- 
though the phase diagrams of these models can be stud- 
ied at weaker couplings (on small lattices away from the 
chiral limit) it would be interesting examine them 
at strong couplings (on large lattices in the chiral limit) . 

Current lattice QCD results at weaker couplings with 
light dynamical quarks suffer from large systematic er- 
rors. Bringing numerical precision, as demonstrated in 
this article, to these studies should be considered an im- 
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portant challenge to pursue in the future. 
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